Paraelectric and ferroelectric order in two-state dipolar fluids 
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Monte Carlo simulations are used to examine cooperative creation of polar state in fluids of two- 
state particles with nonzero dipole in the excited state. With lowering temperature such systems 
undergo a second order transition from nonpolar to polar, paraelectric phase. The transition is 
accompanied by a dielectric anomaly of polarization susceptibility increasing by three orders of 
magnitude. The paraelectric phase is then followed by formation of a nematic ferroelectric which 
further freezes into an fee ferroelectric crystal by a first order transition. A mean-field model of 
phase transitions is discussed. 
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In this Letter we describe Monte Carlo (MC) simula- 
tions of a fluid composed of two-state (TS) soft sphere 
(SS) particles interacting with the dipole-dipole poten- 
tial. This model system, referred to as TS/SS fluid, 
shows a complex phase diagram including the transition 
from a nonpolar to polar, paraelectric phase followed by 
the transition to the ferroelectric phase. The existence 
of ferroelectric order in dipolar systems was first sug- 
gested by Debyei who predicted a transition at y c = 1; 
y — (4ir/9)m 2 p/kT, where m is the dipole moment and 
p is the number density. Onsager«2> and Kirkwood* have 
shifted the transition temperature to zero, y c — > oo. The 
generalized mean-field theory by H0ye and Stell— predicts 
y c = (1 — O) -1 , < O < 1 thus allowing transition at any 
temperature between the Debye and Onsager-Kirkwood 
limits. Recent computer simulations have indicated that 
isotropic ferroelectric phase is possible for polar fluids at 
a non-zero temperature^Sii although the issue of bound- 
ary conditions is still debated^ Most of the discussion of 
spontaneous order in dipolar systems has focused on sys- 
tems with permanent dipoles. Real systems, either of 
molecular or nano-scalc dimension, arc composed of po- 
larizable particles. The present simulations establish the 
existence of spontaneous ferroelectric order in polarizable 
dipolar fluids. 

Many strongly dipolar states are created by in- 
tramolecular separation of charge, either thermal or 
optical, between the the donor and acceptor parts of 
a molecule. Such states, common in chemistry and 
biology, 9 often behave as independent donor-acceptor 
complexes with no cooperativity of charge-transfer tran- 
sitions. A dramatic distinction from this situation is the 
Peierls instability in organic charge-transfer salts where 
charge-separated states are cooperatively created in one- 
dimensional stacks of donor and acceptor units>i& It ap- 
pears that there is no fundamental reason why this type 
of transition should be limited to the crystal phase. We 
show that the cooperative coupling between TS dipoles 
leads to the nonpolar/polar second-order phase transi- 



tion in the liquid phase. 

The Hamiltonian of a fluid of TS/SS particles is 



(uss(jk) - fij 



m fc n fc ) . 



Each individual particle j is characterized by the vacuum 
two-state Hamiltonian Hts(j) with the excited state 
population hj , the vacuum energy gap AI and the mixing 
between the states V. The dipole moment is zero in the 
vacuum ground state and is m in the excited state. The 
excited-state dipoles interact via the dipole-dipole poten- 
tial with Tjk = — VjVfc|rj — r fe | _1 in Eq. (JIJ. Finally, 
the SS repulsion is 



uss{jk) = 4e (a/r jk ) 



12 



(2) 



The model is characterized by the reduced density, p* — 
pa 3 , reduced temperature, T* — kT/e, and reduced 
dipole moment, m* = m/V eer 3 . The reduced parame- 
ters for the TS system are: V* = V/kT and I* = AI/kT. 
The intermolecular interactions are fully characterized by 
two parameters: x — (p*) 4 /T* for SS repulsions^ and y 
for dipolar coupling. 

The exact calculation of the energies of N particles re- 
quires diagonalization of the (2N) 2 matrix. This is still 
prohibitively slow for condensed phase simulations. We 
will therefore use the Hartree decoupling assuming that 
the field acting on a given particle is produced by av- 
erage excited-state populations of other particles^ The 
ground state energy of the fluid E = (JT . ^ g3 \H\ J7 fc ^ gk ) 
is then defined on the wave functions W g j diagonalizing 
the Hamiltonian 



H(j) = HtsU) ~ njirij ■ Kj, 



where 



(3) 
(4) 



is the reaction field of the system dipoles. One gets 
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E = J2 E g (j) + Y,^sW - uddW). (5) 

3 3<k 
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Here, the energy of each particle is 

E g (j) = ±[AI + u Rj -AE j ] (6) 

and the energy of the ground state in vacuum is 
set to be zero of energy. The energy gap AEj — 

[{AI + u Rj ) 2 + AV 2 ] 1/2 in Eq. © is modulated by the 
coupling of the excited-state dipole nx, to the reaction 
field, URj = — m.j ■ Rj. Finally, the excited-state pop- 
ulation 2(rij) = 1 — (AI + ujij)/ AEj defines the dipole 
moment fij = mj(rij) in the dipole-dipole interaction 
potential in Eq. JSJ: UDD(jk) — —fij ■ Tjk ■ Hk- The 
dipole moment of each particle thus fluctuates with the 
instantaneous value of the reaction field. 

The MC simulations of TS/SS fluids were performed 
for the NVT ensemble of N — 256 particles in a cubic 
box with periodic boundary conditions. The long range 
dipolar interactions were accounted for by the reaction 
field method with infinite reaction field dielectric con- 
stant. An MC trial move combines displacement and 
rotation of a particle followed by iterative, self-consistent 
calculation of the fields Rj (i£ Each trial move thus results 
in a new set of dipoles fij. This calculation is the most 
time-consuming part of the simulation protocol. Typical 
simulations were 8x 10 5 cycles long (f 70 h on a single Al- 
pha/740MHz processor), points close to phase transitions 
required simulations with 1.2xl0 6 cycles. 

The appearance of the nonpolar/polar (NP) and po- 
lar/ferroelectric (PF) transitions was monitored by cal- 
culating the polar order parameter P = (fi) /to = 
(Nm) -1 (%2j Mj) an d the usual orientational first-order 
(ferroelectric) and second-order (nematic) parameters, 
S\ and S2, respectively^* The isotropic, nonpolar phase 
is characterized by all order parameters equal to zero. 
The polar state is marked by P ^ 0, whereas P, Si, and 
S2 are all non-zero in the ferroelectric nematic. 
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FIG. 1: Polar order parameter of non-interacting particles 
(line) and of the TS/SS fluid (points) vs the reduced vacuum 
energy gap P = AI/kT; p* = 0.8, m* = 2.0, T* = 0.8. The 
dashed line indicate the point {p)/m = 1/2 at I* =0, the 
dotted line connects the simulation points. 

Spontaneous creation of a non-zero average dipole mo- 
ment is a result of cooperative coupling between induced 
dipoles. The dipole moment (i of an isolated particle 
changes smoothly with the energy gap I* reaching one 
half of its maximum value to at zero gap I* = (Fig. ^| . 



This point corresponds to a purely covalent state of the 
donor-acceptor complex. The ionic character, fi/m ~ 1, 
is achieved at I* < 0, |7*/V*| » 1. In contrast, when the 
particles are allowed to interact in the fluid phase, coop- 
erativity results in a very sharp change of (/i) from zero 
to to at the transition value I^ P > (Fig. [2 Inp = 5-2 
at p* = 0.8, T* = 0.8, to* = 2.0, and V* = -1.0). Ther- 
mal fluctuations of induced dipoles thus reinforce each 
other leading to a negative energy of interaction with the 
reaction field, (ur) < 0, which overrides the positive vac- 
uum energy gap and the polarization energy invested in 
creating the dipole. 

This type of transition was predicted by LoganiSii& 
who suggested that a solution of alkali atoms can show 
a continuous transition from a normal insulator (non- 
polar) to excitonic insulator (paraelectric) phasei 12 i 17 i 18 
Although the atomic dipole is created by hybridizing the 
non-polar atomic states through the transition dipole, 
the basic physics of reaction-field stabilization overriding 
the energy gapi2i2I is common to the present problem 
and the excitonic insulator transition. In is interesting 
to note that non-linear polarizability effects characteris- 
tic of the two-state description might be important for 
the transition to occur since the complex phase diagram 
obtained here for the TS/SS fluid does not exist for the 
fluid of Drude oscillators often used to model atomic and 
molecular polarizability. 21 
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FIG. 2: Order parameters P (squares), Si (circles), and 5*2 
(triangles) of TS/SS fluids vs l/T* (upper panel, m* = 2.0) 
and vs (m*) 2 (lower panel, T* = 1.25) at p* = 0.8, I* = 4.0, 
V* = —1.0. The solid lines show the results of mean-field 
theory calculations [Eqs. (I71- 111H 1 for the order parameter P. 
The Binder parameter^ 1 - (E 4 ) /3{E 2 ) 2 , vs l/T* (upper 
panel) and vs (m*) 2 (lower panel) is shown in the insets. 

The MC simulations shown below consider the varia- 
tion of intermolecular repulsion and attraction through 
T* and to* at constant p* , I*, and V*. Three order pa- 
rameters as functions of T* and to* are shown in Fig. 
The parameter P increases sharply at T^ P /m* NP mark- 
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ing the onset of the polar (paraelectric) phase. Both Si 
and S2 start to increase with decreasing temperature or 
increasing dipole moment resulting in the PF transition 
at Tp F /m* PF . The latter point is defined from the peak 
of dielectric susceptibility (Fig. |3| the simulation data 
corresponding to the change of m* at constant T* are 
not shown). A similar susceptibility peak was observed 
for the PF transition in a fluid of hard sphere Ising spins 
with a square- well attraction^ 

Finally, both S\ and 5*2 have a discontinuous jump at 
Tq = 0.73. The analysis of the density structure factors 
and pair distribution functions indicates that the system 
is in fluid phase above Tq and freezes into the fee fer- 
roelectric at Tq- Further, the phase between Tpp and 
Tc is a ferroelectric nematic as is seen from the analy- 
sis of angular projections of the pair distribution func- 
tion and from the analysis of the position distribution 
functions parallel and perpendicular to the director (not 
shown here) . These functions show no indication of spa- 
tial structure thus ruling out a possibility of smectic or- 
der. No crystallization was achieved by changing m* at 
constant T* in the range of values studied. 
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FIG. 3: Heat capacity C v (upper panel) and polarization sus- 
ceptibility x (middle panel) of the TS/SS fluid and heat ca- 
pacity of the SS fluid (bottom panel), m* = 2.0, I* = 4.0, 
V* = -1.0 for the TS/SS fluid, p* = 0.8 for both the TS/SS 
and SS fluids. MC simulations of the SS fluid were performed 
with N = 2048 particles in the box. The typical simulation 
length was 50 000 cycles. 

The NP transition is characterized by discontinuity of 
the heat capacity C v /Nk = ((SE) 2 )/N(kT) 2 and a very 
substantial (3 orders of magnitude) increase in the dielec- 
tric susceptibility x = p({SM) 2 )/NkT, where <5M is the 
fluctuation of the total dipole of the system of N particles 
(Fig. [3J . At the same time, the total energy is continuous 
pointing to a second-order phase transition (Fig.0J). The 
PF transition is continuous both in energy and heat ca- 



pacity and is manifested by maxima of C v and x (Fig- El- 
Finally, crystallization is a first-order transition with the 
latent heat arising mostly from the change in the dipo- 
lar interaction energy. The continuous character of the 
NP and PF transitions is supported by the value of the 
Binder parameter^ 1 - (E 4 ) /3(E 2 ) 2 , which is very close 
to 2/3 (continuous transition) in the entire range of pa- 
rameters except at crystallization when it takes a dip con- 
sistent with the first order of this transition. The onset of 
the polar phase is also marked by a continuous increase in 
the average reaction filed and a discontinuous jump in its 
variance (Fig. bottom panel). The breakdown of the 
linear response approximation, —kT(uu) = ((Sup) 2 ), at 
the onset of paraelectric phased is characteristic of po- 
larizable systemsiSi 

The temperature of NP transition falls in the region 
of freezing transition of the reference fluid with the re- 
pulsive potential Uss = J2j<k u ss( r jk) which is char- 
acterized by a peak of the heat capacity C v /Nk = 
((SUss) 2 ) /N(kT) 2 (bottom panel in Fig.©. The SS fluid 
crystallizes into an fee lattice below the transition tem- 
perature Tq S as we found from the density structure fac- 
tors in accord with previous reports in the literature i24 
It may be suggested that positional instabilities of the 
reference repulsive potential drive the PF transition of 
the TS/SS fluid. However, dipolar interactions do not 
favor the highly symmetric fee lattice^ and the sys- 
tem stays in the fluid phase in the temperature range 
Tpf < T < Tc crystallizing at Tc- This interpre- 
tation is consistent with the notion advocated by sev- 
eral authors^ that dipoles suppress freezing into a high- 
symmetry lattice due to the anisotropic nature of dipo- 
lar forces. Note that the fee lattice does not necessarily 
represent the ground state since a ferroelectric solid can 
never be strictly cubic. Other stable structures may be 
suppressed in simulations by the cubic shape of the sim- 
ulation cell and periodic boundary conditions. 6 

The energy per particle e = E/NkT depends on the 
fluctuating variable u — un/kT . The mean-field solution, 
neglecting these fluctuations, is given as 

(e) = (P/2){u)-PAe((u)), (7) 

where Ae = AE/kT. The average (u) is affected by the 
macroscopic reaction field 

H« = (47r/3)m/»P5i(l - 0)d (8) 

caused by spontaneous ferroelectric polarization of the 
liquid and the microscopic reaction field R m i C due to 
dipole-dipole correlations. Here, is the mean-field pa- 
rameter of H0ye and Stell^*2& correcting the macroscopic 
field by the account of local dipolar correlations. If one 
assumes that the microscopic correlations are not affected 
by the macroscopic order, (u) becomes 

(u) = 3yP(l - 6)^ - (2/P)u D ((fj)), (9) 

where up>((n)) = —/3(fi) ■ R m i c /2 is the reduced average 
interaction energy between particles carrying the average 
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FIG. 4: Total average energy per particle, E/NkT, the en- 
ergy of dipolar coupling to the reaction field, — (ur), and its 
variance, {(5u R ) 2 ), vs l/T* for the TS/SS fluid: p* = 0.8, 
m* = 2.0, /* =4.0, V* = -1.0. 

dipole (/i) = Pm. Such dipolar interactions are well de- 
scribed by the Pade-truncated perturbation expansion 27 
resulting in the following free energy per particle 

f D kT=-— — (10) 
l + yP 2 I 3 (x)/I 2 {x) 

The perturbation integrals in Eq. I|1(J|I are defined on 
the reference SS fluid as 

h{x) = (27/128 7 r 2 )((p 2 r 6 )- 1 ) ss , 

h(x) = {27/512ir 3 )(p 3 W DDD ) ss , ' ' 



where Wddd is the Axihod- Teller potential^ The av- 
erage over the configuration of the SS fluid, (. . . )sS: nas 
been calculated from MC simulations. hfa) and I3 (x) 
obtained from simulations at different temperatures and 
densities all fall on one universal dependence on x. How- 
ever, the scaling l2(x) oc x~ ' 5 and ^(x) oc x~ ' 75 pre- 
dicted from the statistical average is not seen in the 
simulations. Instead, much weaker power laws, hix) = 
0.1433X" 0138 and I 3 (x) = 0.0314a;- 09915 , fit the data. 
The perturbation integrals obtained at different p* and 
T* are smooth functions of x up to p* = 0.9. 

The paraelcctric order parameter can be calculated 
with the mean-field solution for (u) and the free energy 
obtained by thermodynamic integration of (e) in Eq. (7J . 
The mean-field solution does not reproduce the separa- 
tion of the NP and PF transitions and instead gives a 
direct transition from a non-polar phase to a polar, fer- 
roelectric phase at a single transition point. The calcula- 
tion shown in Fig. [21 is performed for the O parameter in 
Eq. (J5J from the mean-spherical solution for hard dipolar 
spheresA 2 ^ The agreement with simulations for the PF 
transition is very good for both the dependence on tem- 
perature (upper panel in Fig. |2J) and the dependence on 
the dipole moment (lower panel in Fig. 01 . A model in- 
cluding fluctuations of the reaction field 2 *! may result in a 
more accurate description of the NP and PF transitions 
yielding a separate paraelectric phase. 
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